%Model code for 'Expectations, Infections, and Economic Activity' by
%M. Eichenbaum, M. Godinho de Matos, F. Lima, S. Rebelo and M. Trabandt

%Fall 2022. mathias.trabandt@gmail.com

clear all;clc; close all;

load_old_sol=0;

%parameters
par.betta=0.97^(1/52);      %discount factor

par.deltao=1/(13*52);       %prob of dying naturally old (average life expectancy old...)
par.deltay=1/(51*52);       %prob of dying naturally young


par.days=14;                      %average days to either recover or die
par.pidy=0.000409592326139;       %end of sample value of pidy from estimated model %7*0.001/par.days;        %Weekly probability of dying, young
par.piry=7*1/par.days-par.pidy;   %weekly probability of recovering, young
par.pido=0.014335731414868;       %end of sample value of pido from estimated model 7*0.035/par.days;        %Weekly probability of dying, old
par.piro=7*1/par.days-par.pido;   %weekly probability of recovering, old


%beliefs of mortality and recovery rates; only affects FOCs for iy and
%io. If parameters are identical to above ones, no change to results.
do_beliefs=0; %if =1, beliefs of mortality rates; if =0, standard

if do_beliefs==1
    par.pidytilde=0.089052852121211;        %Weekly probability of dying, young
    par.pidotilde=0.428188485888119;        %Weekly probability of dying, old
else
    par.pidytilde=par.pidy;
    par.pidotilde=par.pido;
end
par.pirytilde=7*1/par.days-par.pidytilde;    %Weekly probability of recovering, young
par.pirotilde=7*1/par.days-par.pidotilde;    %Weekly probability of recovering, old


par.phiy=1/(26*1);          %Weeky share of susceptibles that are vaccinated, young
par.phio=1/(26*1);          %Weeky share of susceptibles that are vaccinated, old

par.omy=0.7;                %share of young population, initial
par.omo=1-par.omy;          %share of old pop.

par.pis=1/(1*26);           %reinfection prob.

par.gam=2/3;                %labor share

par.deltak=0.1/52;          %capital depreciation rate

par.rho=1/1.5;              %1/IES

par.z=-1.125;               %constant in utility function, rigged to get VoL of 890k

par.Rnottarget=2.5;                  %R0
par.pi1shrtarget=0.046301363291318;  %share of infections due to consumption
par.pi2shrtarget=0.046301363291318;  %share of infections due to labor

par.inctarget=19000/52;     %income target
par.ntarget=28;             %hours target

%pre-epi steady state
preSy=par.omy;
preSo=1-preSy;
par.nu=1/preSy*(preSo*par.deltao);

preBy=preSy*(par.deltay+par.nu);
preBo=0;

preN=par.ntarget;
preRk=1/par.betta-1+par.deltak;
par.A=(par.inctarget/( ((1-par.gam)/preRk)^((1-par.gam)/par.gam)*preN) )^par.gam;
preK=((1-par.gam)*par.A*preN^par.gam/preRk)^(1/par.gam);
preinc=par.A*preK^(1-par.gam)*preN^par.gam;
preC=preinc-par.deltak*preK;
prew=par.gam*preinc/preN;
preNsy=preN;
preNso=preN;
preCsy=preC;
preCso=preC;
prelamb=preC^(-par.rho);
par.theta=prew*prelamb/preN;
if par.rho==1
    preusy=log(preCsy)-par.theta/2*preNsy^2;
    preuso=log(preCso)-par.theta/2*preNso^2;
else
    preusy=(preCsy^(1-par.rho)-1)/(1-par.rho)-par.theta/2*preNsy^2;
    preuso=(preCso^(1-par.rho)-1)/(1-par.rho)-par.theta/2*preNso^2;
end
prem=preSy*preusy+preSo*preuso;
preU=1/(1-par.betta)*(par.z+prem);
preVoL=preU/prelamb;


%calibrate transmission function params
par.pi1=par.pi1shrtarget*par.Rnottarget*(preSy*par.piry+preSo*par.piro+preSy*par.pidy+preSo*par.pido)/preC^2;
par.pi2=par.pi2shrtarget*par.Rnottarget*(preSy*par.piry+preSo*par.piro+preSy*par.pidy+preSo*par.pido)/preN^2;
par.pi3=par.Rnottarget*(preSy*par.piry+preSo*par.piro+preSy*par.pidy+preSo*par.pido)-par.pi1*preC^2-par.pi2*preN^2;

%solve for endemic steady state
%guesses
By=preBy*100000*1.01;
Bo=0;
N=preN-1;
Ciy=preC-1;
Csy=preC-1;
Nsy=preN-1;
lamtauy=-10;
lamtauo=-10;
lamry=9700/100;
lamro=4600/100;
guess=[By;Bo;N;Ciy;Csy;Nsy;lamtauy;lamtauo;lamry;lamro];


guess=[  71.14
    5.44
    27.92
    272.16
    271.84
    27.97
    -6.91
    -135.32
    102.93
    57.38];

if load_old_sol==1
    load old_sol;
    guess=old_sol;
end


%solve for endemic steady state using fmincon
opts_fmincon=optimoptions('fmincon','Display','iter','ConstraintTolerance',1e-9,'StepTolerance',1e-9,'TolFun',1e-9,...
    'MaxFunctionEvaluations',20000,'UseParallel',false,'FiniteDifferenceStepSize',1e-5,'algorithm','interior-point');

[sol,fval,exitflag]= fmincon(@get_err,guess,[],[],[],[],[],[],@nonlinconstraint,opts_fmincon,par);


if exitflag~=1
    error('fmincon did not solve correctly');
end

old_sol=sol;
save old_sol old_sol;

%get allocations at endemic steady state
[objective,cc,ceq,By,Iy,Io,N,Ciy,Csy,Nsy,lamtauy,lamtauo,Rk,Ry,Ro,Ty,To,Sy,So,K,w,lamb,Cio,Cro,Cry,Niy,Nio,Nry,Nro,C,Nso,Cso,usy,uso,uiy,uio,ury,uro,m,U,dUdSy,dUdSo,dUdRy,dUdRo,dUdIy,dUdIo,lamso,lamsy,lamio,lamiy,lamro,lamry,VoL,Bo]=get_err(sol,par);


%Life expectancy calculations
preLEy=1/(par.deltay+par.nu)*(1+par.nu/par.deltao);
preLEo=1/par.deltao;

%guess LEro;LEio;LEso;LEry;LEiy;LEsy
guess=[13;13;13;27;27;27]*52;
opts_fsolve=optimoptions('fsolve','Display','iter','StepTolerance',1e-9,'TolFun',1e-9);
[sol,fval,exitflag]= fsolve(@get_LE,guess,opts_fsolve,par,To,So,Ty,Sy);
LEro=sol(1);
LEio=sol(2);
LEso=sol(3);
LEry=sol(4);
LEiy=sol(5);
LEsy=sol(6);

%compute percent change in life expectancy
startage_young=39.5; %median age of young at start of epi
startage_old=64.5;
LEatbirth_pre=startage_young+preLEy/52;
LEatbirth_post=startage_young+LEsy/52;
percent_change_LEatbirth=(LEatbirth_post-LEatbirth_pre)/LEatbirth_pre*100;
dlogC_dlogLE=(C/preC-1)*100/percent_change_LEatbirth;




disp(' ');
disp('Steady state results');
disp('--------------------');
disp(' ');
disp('Variable     pre-epi       endemic        percent (or pp) change')
disp('----------------------------------------------------------------');
fprintf('C  (pc)      %6.6f    %6.6f     %6.6f\n',preC,C,(C/preC-1)*100);
fprintf('Cy (pc)      %6.6f    %6.6f     %6.6f\n',preC,(Sy*Csy+Iy*Ciy+Ry*Cry)/(Sy+Iy+Ry),(((Sy*Csy+Iy*Ciy+Ry*Cry)/(Sy+Iy+Ry))/(preC)-1)*100);
fprintf('Co (pc)      %6.6f    %6.6f     %6.6f\n',preC,(So*Cso+Io*Cio+Ro*Cro)/(So+Io+Ro),(((So*Cso+Io*Cio+Ro*Cro)/(So+Io+Ro))/preC-1)*100);
fprintf('N  (pc)      %6.6f     %6.6f      %6.6f\n',preN,N,(N/preN-1)*100);
fprintf('Ny (pc)      %6.6f     %6.6f      %6.6f\n',preN,(Sy*Nsy+Iy*Niy+Ry*Nry)/(Sy+Iy+Ry),(((Sy*Nsy+Iy*Niy+Ry*Nry)/(Sy+Iy+Ry))/preN-1)*100);
fprintf('No (pc)      %6.6f     %6.6f      %6.6f\n\n',preN,(So*Nso+Io*Nio+Ro*Nro)/(So+Io+Ro),(((So*Nso+Io*Nio+Ro*Nro)/(So+Io+Ro))/preN-1)*100);

fprintf('K  (pc)      %6.5f   %6.5f    %6.5f\n',preK,K,((K)/preK-1)*100);
fprintf('w            %6.6f      %6.6f        %6.6f\n',prew,w,(w/prew-1)*100);
fprintf('inc          %6.6f    %6.6f     %6.6f\n',preinc,par.A*K^(1-par.gam)*N^(par.gam),(par.A*K^(1-par.gam)*N^(par.gam)/preinc-1)*100);
fprintf('C/inc        %6.6f      %6.6f        %6.6f\n',preC/preinc,C/(par.A*K^(1-par.gam)*N^(par.gam)),(C/(par.A*K^(1-par.gam)*N^(par.gam))/(preC/preinc)-1)*100);
fprintf('K/inc(ann)   %6.6f      %6.6f       %6.6f\n',preK/preinc/52,K/52/(par.A*K^(1-par.gam)*N^(par.gam)),((K/(par.A*K^(1-par.gam)*N^(par.gam)))/(preK/preinc)-1)*100);
fprintf('lamb         %6.6f      %6.6f       %6.6f\n',prelamb,lamb,(lamb/prelamb-1)*100);
fprintf('Rk           %6.6f      %6.6f        %6.6f\n',preRk,Rk,(Rk/preRk-1)*100);
fprintf('rr(ann)      %6.6f      %6.6f        %6.6f\n',((preRk+1-par.deltak)^52-1)*100,((Rk+1-par.deltak)^52-1)*100,((((Rk+1-par.deltak)^52-1)*100)/(((preRk+1-par.deltak)^52-1)*100)-1)*100);
fprintf('U            %6.5f   %6.5f    %6.5f\n',preU,U,(U/preU-1)*100);
fprintf('m            %6.6f     %6.6f      %6.6f\n',prem,m,(m/prem-1)*100);
fprintf('VoL(mill)    %6.6f      %6.6f        %6.6f\n\n',preVoL/1000000,VoL/1000000,(VoL/preVoL-1)*100);

fprintf('By           %6.6f      %6.6f       %6.6f\n',preBy,By,((By)-preBy));%equals total deaths(natural+Covid)
fprintf('Bo           %6.6f      %6.6f       %6.6f\n\n',preBo,Bo,((Bo)-preBo));%equals total deaths(natural+Covid)


fprintf('Sy           %6.6f      %6.6f       %6.6f\n',preSy,Sy,Sy-preSy);
fprintf('So           %6.6f      %6.6f       %6.6f\n',preSo,So,So-preSo);
fprintf('Iy           %6.6f      %6.6f        %6.6f\n',0,Iy,Iy);
fprintf('Io           %6.6f      %6.6f        %6.6f\n',0,Io,Io);
fprintf('Ry           %6.6f      %6.6f        %6.6f\n',0,Ry,Ry);
fprintf('Ro           %6.6f      %6.6f        %6.6f\n',0,Ro,Ro);
fprintf('Ty           %6.6f      %6.6f        %6.6f\n',0,Ty,Ty);
fprintf('To           %6.6f      %6.6f        %6.6f\n\n',0,To,To);

fprintf('Sy+Iy+Ry     %6.6f      %6.6f        %6.6f\n',preSy,Sy+Iy+Ry,Sy+Iy+Ry-preSy);
fprintf('So+Io+Ro     %6.6f      %6.6f        %6.6f\n\n',preSo,So+Io+Ro,So+Io+Ro-preSo);


fprintf('Dweek        %6.6f      %6.6f        %6.6f\n',0,(par.pidy*Iy+par.pido*Io),(par.pidy*Iy+par.pido*Io));
fprintf('Dann         %6.6f      %6.6f        %6.6f\n\n',0,(par.pidy*Iy+par.pido*Io)*52,(par.pidy*Iy+par.pido*Io)*52);


fprintf('Csy          %6.6f    %6.6f     %6.6f\n',preC,Csy,(Csy/preC-1)*100);
fprintf('Cso          %6.6f    %6.6f     %6.6f\n',preC,Cso,(Cso/preC-1)*100);
fprintf('Ciy          %6.6f    %6.6f      %6.6f\n',preC,Ciy,(Ciy/preC-1)*100);
fprintf('Cio          %6.6f    %6.6f      %6.6f\n',preC,Cio,(Cio/preC-1)*100);
fprintf('Cry          %6.6f    %6.6f      %6.6f\n',preC,Cry,(Cry/preC-1)*100);
fprintf('Cro          %6.6f    %6.6f      %6.6f\n\n',preC,Cro,(Cro/preC-1)*100);

fprintf('Nsy          %6.6f     %6.6f      %6.6f\n',preN,Nsy,(Nsy/preN-1)*100);
fprintf('Nso          %6.6f     %6.6f      %6.6f\n',preN,Nso,(Nso/preN-1)*100);
fprintf('Niy          %6.6f     %6.6f      %6.6f\n',preN,Niy,(Niy/preN-1)*100);
fprintf('Nio          %6.6f     %6.6f      %6.6f\n',preN,Nio,(Nio/preN-1)*100);
fprintf('Nry          %6.6f     %6.6f      %6.6f\n',preN,Nry,(Nry/preN-1)*100);
fprintf('Nro          %6.6f     %6.6f      %6.6f\n\n',preN,Nro,(Nro/preN-1)*100);


fprintf('usy          %6.6f     %6.6f       %6.6f\n',prem,usy,(usy/prem-1)*100);
fprintf('uso          %6.6f     %6.6f      %6.6f\n',prem,uso,(uso/prem-1)*100);
fprintf('uiy          %6.6f     %6.6f       %6.6f\n',prem,uiy,(uiy/prem-1)*100);
fprintf('uio          %6.6f     %6.6f       %6.6f\n',prem,uio,(uio/prem-1)*100);
fprintf('ury          %6.6f     %6.6f       %6.6f\n',prem,ury,(ury/prem-1)*100);
fprintf('uro          %6.6f     %6.6f       %6.6f\n\n',prem,uro,(uro/prem-1)*100);


fprintf('lamtauy      %6.6f        %6.6f     %6.6f\n',NaN,lamtauy,NaN);
fprintf('lamtauo      %6.6f        %6.6f   %6.6f\n',NaN,lamtauo,NaN);
fprintf('lamsy        %6.6f        %6.6f  %6.6f\n',NaN,lamsy,NaN);
fprintf('lamso        %6.6f        %6.6f   %6.6f\n',NaN,lamso,NaN);
fprintf('lamiy        %6.6f        %6.6f  %6.6f\n',NaN,lamiy,NaN);
fprintf('lamio        %6.6f        %6.6f   %6.6f\n',NaN,lamio,NaN);
fprintf('lamry        %6.6f        %6.6f  %6.6f\n',NaN,lamry,NaN);
fprintf('lamro        %6.6f        %6.6f   %6.6f\n\n',NaN,lamro,NaN);

fprintf('LEro (years)  %6.6f      %6.6f   %6.6f\n',preLEo/52,LEro/52,NaN);
fprintf('LEio (years)  %6.6f      %6.6f   %6.6f\n',preLEo/52,LEio/52,NaN);
fprintf('LEso (years)  %6.6f      %6.6f   %6.6f\n',preLEo/52,LEso/52,NaN);
fprintf('LEry (years)  %6.6f      %6.6f   %6.6f\n',preLEy/52,LEry/52,NaN);
fprintf('LEiy (years)  %6.6f      %6.6f   %6.6f\n',preLEy/52,LEiy/52,NaN);
fprintf('LEsy (years)  %6.6f      %6.6f   %6.6f\n\n',preLEy/52,LEsy/52,NaN);

fprintf('LE@birth(yrs) %6.6f      %6.6f    %6.2f (prct change) \n',LEatbirth_pre,LEatbirth_post,percent_change_LEatbirth);
fprintf('dlogC_dlogLE  %6.6f         %6.6f        %6.3f \n\n',NaN,NaN,dlogC_dlogLE);
 

 
disp('Parameters calibrated in pre-epi steady state (kept fixed)');
fprintf('A              %6.6f       %6.6f      %6.6f\n',par.A,par.A,0);
fprintf('theta          %6.6f       %6.6f      %6.6f\n',par.theta,par.theta,0);
fprintf('nu             %6.6f       %6.6f      %6.6f\n\n',par.nu,par.nu,0);


fprintf('Age young (data, years)  %6.6f   \n',startage_young);
fprintf('Age old (data, years)    %6.6f   \n\n',startage_old);


function [objective,cc,ceq,By,Iy,Io,N,Ciy,Csy,Nsy,lamtauy,lamtauo,Rk,Ry,Ro,Ty,To,Sy,So,K,w,lamb,Cio,Cro,Cry,Niy,Nio,Nry,Nro,C,Nso,Cso,usy,uso,uiy,uio,ury,uro,m,U,dUdSy,dUdSo,dUdRy,dUdRo,dUdIy,dUdIo,lamso,lamsy,lamio,lamiy,lamro,lamry,VoL,Bo]=get_err(guess,par)


By=abs(guess(1))/100000;
Bo=(guess(2))/100000;
N=abs(guess(3));
Ciy=abs(guess(4));
Csy=abs(guess(5));
Nsy=abs(guess(6));
lamtauy=guess(7);
lamtauo=guess(8);
lamry=abs(guess(9))*100;
lamro=abs(guess(10))*100;


Rk=1/par.betta-1+par.deltak;

Iy=(By-(par.nu+par.deltay)*par.omy)/((1-par.nu-par.deltay)*par.pidy);

Ty=(1-(1-par.piry-par.pidy)*(1-par.deltay-par.nu))/(1-par.deltay-par.nu)*Iy;

Ry=(1-par.nu-par.deltay)/(par.nu+par.deltay+par.pis+par.phiy*(1-par.nu-par.deltay))*(Iy*par.piry+par.phiy*(par.omy-Iy));

Sy=par.omy-Iy-Ry;


Io=(Bo+par.nu*(Sy+(1-par.pidy)*Iy+Ry)-par.omo*par.deltao)/((1-par.deltao)*par.pido);

To=((1-(1-par.piro-par.pido)*(1-par.deltao))*Io-par.nu*(Ty+Iy*(1-par.piry-par.pidy)))/(1-par.deltao);

Ro=(par.nu*(Ry+Iy*par.piry+par.phiy*Sy)+(1-par.deltao)*(Io*par.piro+par.phio*(par.omo-Io)))/(par.deltao+par.pis+par.phio*(1-par.deltao));

So=par.omo-Ro-Io;


K=(((1-par.gam)*par.A*N^par.gam)/Rk)^(1/par.gam);
w=par.gam*par.A*K^(1-par.gam)*N^(par.gam-1);

lamb=Ciy^(-par.rho);

Cio=Ciy;
Cro=Ciy;
Cry=Ciy;
Niy=w*lamb/par.theta;
Nio=Niy;
Nry=Niy;
Nro=Niy;
C=par.A*K^(1-par.gam)*N^par.gam-par.deltak*K;
Nso=1/So*(N-Iy*Niy-Io*Nio-Ry*Nry-Ro*Nro-Sy*Nsy);
Cso=1/So*(C-Iy*Ciy-Io*Cio-Ry*Cry-Ro*Cro-Sy*Csy);

if par.rho==1
    usy=log(Csy)-par.theta/2*Nsy^2;
    uso=log(Cso)-par.theta/2*Nso^2;
    
    uiy=log(Ciy)-par.theta/2*Niy^2;
    uio=log(Cio)-par.theta/2*Nio^2;
    
    ury=log(Cry)-par.theta/2*Nry^2;
    uro=log(Cro)-par.theta/2*Nro^2;
else
    usy=(Csy^(1-par.rho)-1)/(1-par.rho)-par.theta/2*Nsy^2;
    uso=(Cso^(1-par.rho)-1)/(1-par.rho)-par.theta/2*Nso^2;
    
    uiy=(Ciy^(1-par.rho)-1)/(1-par.rho)-par.theta/2*Niy^2;
    uio=(Cio^(1-par.rho)-1)/(1-par.rho)-par.theta/2*Nio^2;
    
    ury=(Cry^(1-par.rho)-1)/(1-par.rho)-par.theta/2*Nry^2;
    uro=(Cro^(1-par.rho)-1)/(1-par.rho)-par.theta/2*Nro^2;
end

m=Sy*usy+So*uso+Iy*uiy+Io*uio+Ry*ury+Ro*uro;

U=1/(1-par.betta)*(par.z+m);

dUdSy=usy;
dUdSo=uso;

dUdIy=uiy;
dUdIo=uio;

dUdRy=ury;
dUdRo=uro;

lamso=1/(1-par.betta*(1-par.deltao)*(1-par.phio))*(par.betta*dUdSo+par.betta*lamb*(w*Nso-Cso)+par.betta*lamro*par.phio*(1-par.deltao)+par.betta*lamtauo*(1-par.phio)*(par.pi1*Cso*(Iy*Ciy+Io*Cio)+par.pi2*Nso*(Iy*Niy+Io*Nio)+par.pi3*(Iy+Io)));
lamsy=1/(1-par.betta*(1-par.deltay-par.nu)*(1-par.phiy))*(par.betta*dUdSy+par.betta*lamb*(w*Nsy-Csy)+par.betta*lamso*par.nu*(1-par.phiy)+par.betta*lamry*par.phiy*(1-par.deltay-par.nu)+par.betta*lamro*par.nu*par.phiy+par.betta*lamtauy*(1-par.phiy)*(par.pi1*Csy*(Iy*Ciy+Io*Cio)+par.pi2*Nsy*(Iy*Niy+Io*Nio)+par.pi3*(Iy+Io)));

lamio=1/(1-par.deltao)*(lamso*(1-par.deltao)+lamtauo);
lamiy=1/(1-par.deltay-par.nu)*(lamsy*(1-par.deltay-par.nu)+lamso*par.nu+lamtauy-lamio*par.nu);

VoL=U/lamb;


err=NaN*zeros(10,1);
err(1)=(1-par.phiy)*(par.pi1*Sy*Csy*(Iy*Ciy+Io*Cio)+par.pi2*Sy*Nsy*(Iy*Niy+Io*Nio)+par.pi3*Sy*(Iy+Io))-Ty;
err(2)=(1-par.phio)*(par.pi1*So*Cso*(Iy*Ciy+Io*Cio)+par.pi2*So*Nso*(Iy*Niy+Io*Nio)+par.pi3*So*(Iy+Io))-To;
err(3)=Csy^(-par.rho)-lamb+lamtauy*par.pi1*(1-par.phiy)*(Iy*Ciy+Io*Cio);
err(4)=Cso^(-par.rho)-lamb+lamtauo*par.pi1*(1-par.phio)*(Iy*Ciy+Io*Cio);
err(5)=-par.theta*Nsy+w*lamb+lamtauy*par.pi2*(1-par.phiy)*(Iy*Niy+Io*Nio);
err(6)=-par.theta*Nso+w*lamb+lamtauo*par.pi2*(1-par.phio)*(Iy*Niy+Io*Nio);

err(7)=par.betta*dUdIy+par.betta*lamb*(w*Niy-Ciy)-lamiy+par.betta*lamiy*(1-par.pirytilde-par.pidytilde)*(1-par.deltay-par.nu)+par.betta*lamio*(1-par.pirytilde-par.pidytilde)*par.nu+par.betta*lamry*par.pirytilde*(1-par.deltay-par.nu)+par.betta*lamro*par.pirytilde*par.nu;
err(8)=par.betta*dUdIo+par.betta*lamb*(w*Nio-Cio)-lamio+par.betta*lamio*(1-par.pirotilde-par.pidotilde)*(1-par.deltao)+par.betta*lamro*par.pirotilde*(1-par.deltao);

err(9)=-lamro+1/(1-par.betta*(1-par.deltao-par.pis))*(par.betta*dUdRo+par.betta*lamb*(w*Nro-Cro)+par.betta*lamso*par.pis);
err(10)=-lamry+1/(1-par.betta*(1-par.deltay-par.nu-par.pis))*(par.betta*dUdRy+par.betta*lamb*(w*Nry-Cry)+par.betta*lamsy*par.pis+par.betta*lamro*par.nu);


%ceq=0, equality constraint
ceq=err;


%%cc<=0, inequality constraint.
%cc=NaN*zeros(2,1);
cc(1)=-Iy;
cc(2)=-Io;

cc(3)=-Sy;
cc(4)=-So;

cc(5)=-Ry;
cc(6)=-Ro;

cc(7)=-Ty;
cc(8)=-To;

cc(9)=Iy-par.omy;
cc(10)=Io-par.omo;


cc(11)=Sy-par.omy;
cc(12)=So-par.omo;

cc(13)=Ry-par.omy;
cc(14)=Ro-par.omo;

cc(15)=Ty-par.omy;
cc(16)=To-par.omo;

cc(17)=-(By-(par.nu+par.deltay)*par.omy);
%cc=[];



%objective to be minimized (trivial)
objective=0;


end





function [cc,ceq] = nonlinconstraint(guess,par)
%get value of nonlinear constraint   must hold
%with equality (ceq)
[~,cc,ceq]=get_err(guess,par);
end




%function to calculate LE's
function [chk]=get_LE(guess,par,To,So,Ty,Sy)

LEro=abs(guess(1));
LEio=abs(guess(2));
LEso=abs(guess(3));
LEry=abs(guess(4));
LEiy=abs(guess(5));
LEsy=abs(guess(6));


chk(1)=-LEro+par.deltao+(1-par.deltao-par.pis)*(1+LEro)+par.pis*(1+LEso);
chk(2)=-LEio+par.pido+par.piro*par.deltao+(1-par.piro-par.pido)*par.deltao+par.piro*(1-par.deltao)*(1+LEro)+(1-par.piro-par.pido)*(1-par.deltao)*(1+LEio);
chk(3)=-LEso+(1-To/So-par.phio)*(1-par.deltao)*(1+LEso)+To/So*(1-par.deltao)*(1+LEio)+par.phio*(1-par.deltao)*(1+LEro)+par.deltao;

chk(4)=-LEry+par.deltay+par.nu*(1+LEro)+par.pis*(1+LEsy)+(1-par.deltay-par.nu-par.pis)*(1+LEry);
chk(5)=-LEiy+par.pidy+par.piry*par.deltay+(1-par.piry-par.pidy)*par.deltay+par.piry*(1-par.deltay-par.nu)*(1+LEry)+par.piry*par.nu*(1+LEro)+(1-par.piry-par.pidy)*(1-par.deltay-par.nu)*(1+LEiy)+(1-par.piry-par.pidy)*par.nu*(1+LEio);
chk(6)=-LEsy+Ty/Sy*(1-par.deltay-par.nu)*(1+LEiy)+Ty/Sy*par.nu*(1+LEio)+(1-Ty/Sy-par.phiy)*(1-par.deltay-par.nu)*(1+LEsy)+(1-Ty/Sy-par.phiy)*par.nu*(1+LEso)+par.phiy*(1-par.deltay-par.nu)*(1+LEry)+par.phiy*par.nu*(1+LEro)+par.deltay;

end

